#!/usr/bin/perl
 use strict;
# use warnings;
use Cwd;
use lib getcwd."/bin/bioperl-1.5.2_102"; # add local::
use lib getcwd."/lib";
use Bio::TreeIO;
##################################################### INI File
# my ($USAGE) =  "\n\n\t****** USAGE of preprocess.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> \n\n\n\n";
# die "$USAGE" if ((scalar(@ARGV))< 1);
# my $option_file = $ARGV[0];
my $option_file = "panam.ini";
chomp($option_file);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;
my $path = "bin/uclust3.0.617";
my $us_version = "3.0.617";
####
my $lib = "lib";
my $fast = "$lib/FastTree";
#my $align = "$lib/hmmalign";
my $Reference = "Reference" ;
my $build = "$lib/hmmbuild";
####
my $dom;
my %dom;
my $aff_norm;


while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


#		if($line=~m/QUERY_SEQUENCES/gi) {
#			$line=~s/QUERY_SEQUENCES//gi;
#			$line=~s/\t//gi;
#			$line=~s/^\s+//gi;
#			$line=~s/\s+$//gi; 
#			if($line) {
#				$query_seq = $line; 
#			}
#			chomp ($query_seq);
#			# See if the query sequences file exist and contains something
#			die "\n\n\tCannot find query sequence file: $query_seq. Check panam.ini.\n\n" unless (-e $query_seq);
#			die "\n\n\tQuery sequence file, $query_seq, appears to be empty!\n\n" if (-z $query_seq);
#			
#			my $ligne1= `sed -n '1p' $query_seq`;
#			my $ligne2= `sed -n '2p' $query_seq`;	
#			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[atcguATCGU]/)) {
#				die "\n\n\tYour file $query_seq does not seem to contain sequences in fasta format.\n\n";
#			}		
#		}

		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}

#		if($line=~m/REFERENCE_BASE/gi) {
#			$line=~s/REFERENCE_BASE//gi;
#			$line=~s/\t//gi;
#			$line=~s/^\s+//gi;
#			$line=~s/\s+$//gi;
#			if($line) {
#				$Reference = $line; 
#			}
#		}

		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}

		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;				
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


print "$Reference\n";

########################### ajout position

my $forw_regexp; my $reve_regexp;
my @tab_forw=split(//,$seq_F);
$forw_regexp=regexp(@tab_forw);

my $rev=reverse($seq_R);
my @tab_reve=split(//,$rev);
$reve_regexp=regexp(@tab_reve);
$reve_regexp=~tr/ACGURMBD/UGCAYKVH/;

sub regexp{
	my $regexp;
	foreach (@_){
		if ($_ eq "T"){
			$regexp.="U";
		}
		elsif ($_ eq "R"){
			$regexp.="[AGR]";
		}
		elsif ($_ eq "Y"){
			$regexp.="[CUY]";
		}
		elsif ($_ eq "M"){
			$regexp.="[CAM]";
		}
		elsif ($_ eq "K"){
			$regexp.="[UGK]";
		}
		elsif ($_ eq "W"){
			$regexp.="[UAW]";
		}
		elsif ($_ eq "S"){
			$regexp.="[CGS]";
		}
		elsif ($_ eq "B"){
			$regexp.="[CUGB]";
		}
		elsif ($_ eq "D"){
			$regexp.="[AUGD]";
		}
		elsif ($_ eq "H"){
			$regexp.="[AUCH]";
		}
		elsif ($_ eq "V"){
			$regexp.="[ACGV]";
		}
		elsif ($_ eq "N"){
			$regexp.="[ACGUN]";
		}
		else{
			$regexp.=$_;
		}
	}
	return $regexp;
}

############################################# trimming profiles ######################################################

if (!(defined $primF) and (defined $primR))  {
	die "\n\n\t Primer forward not defined. Check panam.ini.\n\n";
}

if (defined $primF) {
	die "\n\n\t Primer reverse not defined. Check panam.ini.\n\n" unless (defined $primR);
	die "\n\n\t Primers' sequences missed. Check panam.ini.\n\n" unless ((defined $seq_F) and (defined $seq_R));
	die "\n\n\t Primers' sequences seem to be in the wrong format. Check panam.ini.\n\n" unless (($seq_F =~ /[ATCGUatcgu]/) and ($seq_R =~ /[ATCGUatcgu]/));
	$trim = 1;
}

my @tax=();
my $Profiles;
if ($trim == 1) {
	$Profiles = "Reference/Profiles_".$primF."_".$primR;
	$Profiles = "Reference/Trimmed_Profiles";
	if ((-d $Profiles)) {
		print "Trimmed profiles exist in $Profiles/\n";

	}
	
	else {
		print "\n\tTrimming profiles ...\n"; 

		###############"""récuperation des positions
		my %list_profils;
		my $acces; 
		my $taxo;
# 		my @tax = ();
		my $rep = "Reference/Profiles/fasta/";
		my @listefic = <$rep*.fasta>;
			
		my $pos_forward; my $pos_reverse;
			
		foreach my $nn (@listefic) { 
			foreach my $domain (keys %dom){
				if ($nn =~ /$domain/) { 
					if ($nn=~ /Reference\/Profiles\/fasta\/(.*?).fasta/) { 
						$taxo = $1 ; 
						push (@tax, $taxo);
					}	
					# debug
					# print $taxo."\n";
					open (my $profils,"<Reference/Profiles/fasta/$taxo.fasta") || die "can not open file";
					while (<$profils>){
						chomp $_;
						if ($_=~/^>(.*)$/){
							$acces=$1;
						}
						else{
							$_=~s/ //g;
							$list_profils{$taxo}{$acces}.=$_;
						}
					}
				close $profils;
				}
			}
		}
	
		#print "Taxo : ".join("\n", @tax)."\n\n";
			
		my $pos_reverse_som=0; my $pos_forward_som=0; my $if =0; my $ir=0; 
		foreach my $taxo (@tax){
			my $pos_forward;my $pos_reverse;
			foreach my $acces (keys %{$list_profils{$taxo}}){
				my $sequence=$list_profils{$taxo}{$acces};
		
				$sequence=~s/-//g;
				$sequence=~s/\.//g;
		
				if ($sequence=~/$forw_regexp/){
					$if++;
					$pos_forward=length($`) ;
					$pos_forward_som+=$pos_forward;
				}
				if ($sequence=~/$reve_regexp/){
					$ir++;
					$pos_reverse=length($`) + length($&);
					$pos_reverse_som+= $pos_reverse;
				}
						
				if (defined($pos_reverse) && defined($pos_forward)){
					last;	
				}
			}
		}

		if ($if==0) { die("\#\#ABORT\#\#FORWARD_PRIMER_NOT_FOUND"); }
		if ($ir==0) { die("\#\#ABORT\#\#REVERSE_PRIMER_NOT_FOUND"); }
					
		my $posF = int ($pos_forward_som/$if); 
		my $posR = int ($pos_reverse_som/$ir);

		#############fin recuperation des positions

		trim ($primF, $primR, $posF, $posR);
		print "\n\tTrimmed profiles have been generated in ".$Profiles."\n\n";
	}
}
############################# Trim subroutine
sub trim {
	my $posF; my $posR; my $primF; my $primR; my $fasta_file; my $hmm_file;
	($primF, $primR, $posF, $posR) = @_;

	my $nom = $primF."_".$primR ;
	
	`mkdir Reference/Trimmed_Profiles/`;
	`mkdir Reference/Trimmed_Profiles/hmmprofil/`;
	`mkdir Reference/Trimmed_Profiles/fasta/`;
	
	open (C, "Reference/seq_model_trim") || die "cannot open file";
	my %corresp;
	while (<C>) {
		my @tab = split (/\t/, $_);
		chomp ($tab[0]); chomp ($tab[1]);
		$corresp{$tab[0]} = $tab[1];
	}
	close C;
	
	foreach my $n (@tax) {
		open (F, "Reference/Profiles/fasta/".$n.".fasta") || die "can not open file";
		my %seq ; my $a;
		while (<F>) {
			if ($_=~ />(.*?)\s+/ ) {
				$a = $1;
			}
			else {
				chomp($_);
				$seq{$a}.= $_;
			}
		}
		close F;
		
		my $seq = $corresp{$n};
		
		my @tab = split (//, $seq{$seq}) ;
		
		my $pos_reel =1;
		my $pos_nuc =1;
		my $pos_nuc_fin;
		
		my $posF_gaped;
		my $posR_gaped;
		
		foreach my $e (@tab) {
		
			if (($e ne "-") and ($e ne ".")) {
				$pos_nuc ++;
			}
			if ($pos_nuc == $posF) {
				$posF_gaped = $pos_reel;	
			}
			
			elsif ($pos_nuc == $posR) {
				$posR_gaped = $pos_reel;	
				last;
			}
	
			$pos_reel++;
		}
			
		if (!(defined $posR_gaped)) {$posR_gaped = $pos_reel}
			
		my $length = $posR_gaped - $posF_gaped;
		
		#open (A, ">Reference/Profiles_".$nom."/fasta/".$n.".fasta") || die "can not open file"; 
		open (A, ">$Reference/Trimmed_Profiles/fasta/".$n.".fasta") || die "can not open file"; 
		foreach my $k (keys %seq) {
			my $sub = substr ($seq{$k}, $posF_gaped, $length);
			if ($sub=~ /[ATCGUatcgu]/) {
				print A ">$k\n$sub\n"
			}
		}
		
		close A;
			
		#$fasta_file = "Reference/Profiles_".$nom."/fasta/".$n.".fasta" ;
		#`bin/./hmmer-2.3.2/src/hmmbuild -g "Reference/Profiles_"$nom"/hmmprofil/$n".hmmprofil $fasta_file` ;
		$fasta_file = "$Reference/Trimmed_Profiles/fasta/".$n.".fasta" ;
		`lib/hmmbuild -g "$Reference/Trimmed_Profiles/hmmprofil/$n".hmmprofil $fasta_file` ;
			
		#$hmm_file = "Reference/Profiles_".$nom."/hmmprofil/".$n.".hmmprofil";
		$hmm_file = "$Reference/Trimmed_Profiles/hmmprofil/".$n.".hmmprofil";
	}
	return ($fasta_file , $hmm_file);
}

 
if ($trim == 0) {
	$Profiles= "Reference/Profiles";
}
























